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Abstract . Computer simulations based on the molecular dynamics (MD) technique have been used to study the 
mechanisms and kinetics of free evaporation from crystalline and molten forsterite (i. e., Mg2Si(>4) on an atomic 
level. The interatomic potential employed for these simulations reproduces the energetics of bonding in forsterite 
and in gas-phase MgO and Si02 reasonably accurately. Results of the simulation include predicted evaporation rates, 
diffusion rates, and reaction mechanisms for Mg2Si04(s or 1) 2Mg(g) + 20(g) + SiCtyg). 


Introduction . Hashimoto [1] showed that evaporation and condensation kinetics played a role in the cosmochemistry 
of the early solar nebula. Davis et al. [2] demonstrated based on the results of evaporation experiments on forsterite 
crystals and melt that the so-called FUN CAI's (fractionated and unknown nuclear effect; calcium-, aluminum-rich 
refractory inclusions) in carbonaceous chondrites could have obtained anomalous isotopic signatures from kinetic 
isotope effects during non-equilibrium evaporation. In their experiments, however, isotope fractionations were 
significant only for molten Mg2SiC>4 and not crystalline forsterite. They explained this as due to the fact that the 
rate of evaporation from the solid is much faster than the rate of diffusion, so no isotopic fractionation is possible. 
In contrast, diffusion and evaporation rates operate on similar time scales in the molten phase, so a significant vapor- 
liquid isotopic fractionation occurs because light isotopes from the interior of the melt diffuse to the surface as fast 
as they are evaporated. This study was conducted to investigate possible mechanisms of evaporation from forsterite 
crystals and melts. Such simulations can lead to insights into the main evaporating species (i.e., SiO(g) vs. 
Si02(g)) and into the factors influencing the magnitude of kinetic isotope effects accompanying evaporation. 



Methodology . MD simulations calculate the motion of particles in an atomic system with Newton s equations of 
motion employed over small time increments, usually 10'^ s. The extremely small time steps are required to 
conserve the energy of the system but this causes the real time duration of the simulation to be short (i. e., l(H^ to 
10’ seconds). Structural, thermodynamic, and kinetic properties of the system can be calculated with the Ml) 
technique. (For a review, see [3, 4]). flic MD system in this study consisted of 1680 atoms (240 Mg 2 Si (>4 
formula units) with initial crystal positions determined from the experimental crystal coordinates. Mg2Si<>4 melts 
were generated from the simulations of (51 with 1280 atoms in the MD cell. Two-dimensional periodic boundary 
conditions were employed along the a- and b-axes, approximating an infinite layer with a thickness along the c-axis 
of 24 A in the crystalline phase (Fig. 1) and 16 A in the molten phase (Fi^. 2). The potential of Leinenweber and 
Navrotsky [6] was used to model the interactions of Si 4+ , Mg^+, and ()*•“ ions within the evaporating forsterite. 
We chose this potential because it accurately reproduces mineral structures in the Mg - Si - () system |6], and also 
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does an accurate job of fitting the energetics of neutral gas phase molecules [7]. Evaporation rates were determined 
from simulations at 3000, 4000 and 5000 K in the crystalline phase and 4000, 5000, and 6000 K in the molten 
phase. The positions of all particles are known as a function of time, so the evaporation rate is the number of ions 
leaving the mineral surface per unit time per unit surface area. Atomic masses are read in as part of the input data 
set; hence, it is possible to isotopically "dope" the simulated mineral with ^O, ^Mg, and ™Si. Relative rates of 
evaporation versus diffusion (calculated with the Einstein equation) can be compared in the crystalline and melt 
phases. 

Results and Discussion . Due to the short durations of the MD simulations, the evaporation and diffusion rates are 
most accurately determined at temperatures much higher than normally employed in the laboratory. However, the 
short time duration also allows crystalline systems to be studied at temperatures far above the melting point because 
the system remains crystalline metastably [8]. Evaporation rates predicted from the simulations tire on the order of 1 
to 10 mol/cm^/sec at 3000 to 6000 K in both crystalline and molten Mg 2 Si 04 similar to values obtained from 
extrapolation of the experimental data [1]. Evaporation from crystals occurs more rapidly than from melts at the 
same temperature in these simulations, but diffusion is more rapid in the melt phase. This is also evidenced in the 
melt simulation by mixing of isotopes initially placed in separate regions of the cell (Fig. 2): whereas, in the 
simulations of crystalline Mg 2 SiC> 4 , the isotopes do not mix significantly from one crystal plane to the next (Fig. 
1). These results are consistent with Hashimoto's [1] interpretation of the experimental data. 


as heavy isotopes ( l8 (), 26 Mg, and ? 0 Si), initially in the right-hand side of the cell, mix with light 
isotopes on the left. (Key in Fig. 1). 

The reaction mechanism predicted in the MD simulations is somewhat different from the experimental 
interpretation, however. Although the ionic nature of the interatomic potential may be somewhat inaccurate for 
modelling neutral atoms, The species leaving the surface are dominated by Mg^ + and [SiOq] 4 *' ions rather than Mg 
and O atoms and Si (>2 molecules . Such speciation may contribute to the larger isotopic fractionation of Mg during 
evaporation of Mg 2 SiC >4 melts compared to O and Si [2j because the larger [SiOql^ - anion will be less affected by 
small mass differences due to isotopic substitutions comp;ired to the atoms and molecular Si< >2 *- Gas-phase 
reactions could subsequently lead to conversion of 2Mg 2+ + [ Si< > 4 — > Si( >2 + 2< ) + 2Mg. 
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